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1 Abstract 

Population synthesis is used to model the number of neutron stars in glob- 
ular clusters that are observed as low-mass X-ray sources and millisecond 
radio pulsars. The dynamical interactions between binary and single stars 
in a cluster are assumed to take place only with a continuously replenished 
"background" of single stars whose properties keep track of the variations in 
parameters of the cluster as a whole and the evolution of single stars. We 
use the hypothesis that the neutron stars forming in binary systems from 
components with initial masses of ~ 8 — 12M during the collapse of degen- 
erate O-Ne-Mg cores through electron captures do not acquire a high space 
velocity. The remaining neutron stars (from single stars with masses > 8M 
or from binary components with masses > 12M Q ) are assumed to be born 
with high space velocities. According to this hypothesis, a sizeable fraction 
of the forming neutron stars remain in globular clusters (about 1000 stars in 
a cluster with a mass of 5 • 1O 5 M ). The number of millisecond radio pulsars 
forming in such a cluster in the case of accretion- driven spinup in binary sys- 
tems is found to be 10, in agreement with observations. Our modeling also 
reproduces the observed shape of the X-ray luminosity function for accreting 
neutron stars in binary systems with normal and degenerate components and 
the distribution of spin periods for millisecond pulsars. 

Key words: binary X-ray sources, neutron stars, millisecond pulsars, globular 
clusters. 



2 INTRODUCTION 



Neutron stars (NSs) in globular clusters (GCs) are observed by currently 
available X-ray astronomical methods as accreting compact objects in low- 
mass X-ray binaries (including X-ray bursters) and by radio-astronomical 
methods as millisecond radio pulsars that are generally members of binary 
systems. There is an evolutionary connection between the two classes of 
objects - millisecond radio pulsars are accretion-driven spin-up products of 
old NSs in close low-mass X-ray binaries (this idea was first put forward 
by Bisnovatyi-Kogan and Romberg (1974); see also the review article by 
Bisnovatyi- Kogan (2006); for the most recent observational confirmations, 
see, e.g., Bogdanov et al. (2005)). The standard scenario for the evolution 
of low- mass Xray binaries and millisecond pulsars (see, e.g., Bhattacharya 
and van den Heuvel 1991) in dense GCs is complemented by the effect of 
dynamical interactions between stars - new binaries formed through tidal in- 
teractions, captures, component exchanges, etc. are added to the originally 
existing binary systems. 

Observational data clearly show an increased concentration of X-ray sources 
and binary millisecond pulsars in GCs (see, e.g., Heinke et al. 2005; Grindlay 
2004; Lorimer 2005), thereby confirming that dynamical effects are crucially 
important for describing their formation and evolution (Pooley et al.2003). 

NSs are formed through the collapse of the cores of stars with main-sequence 
masses above a certain value 8-10 M Q , determined by the chemical com- 
position and, possibly, other parameters, that finished their thermonuclear 
evolution. Depending on the evolutionary status of the star at the time of 
its collapse, the NS birth through core collapse is accompanied by a type-II 
or type-Ibc supernova explosion and by the ejection of a significant presu- 
pernova mass in the form of a shell. Data on young NSs observed as radio 
pulsars are also indicative of their high space velocities with amean value of 
~ 400 km s _1 (see, e.g., Lorimer (2005) and references therein; Hobbs et al. 
2005). The hypothesis about a kick velocity acquired by a NS during collapse 
is commonly invoked to account for such high space velocities. The physical 
mechanism underlying this velocity is unknown and there are many specu- 
lations and models on this score; many of them were listed and discussed in 
the review by Lai (2001). Note that deducing the shape of the distribution 
of kicks from observations is coupled with difficult allowance for selection 
effects (Tutukov et al. 1984) and it is believed (see, e.g., Ib en and Tutukov 
1996) that nature does without kicks and all of the high space velocities of 
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radio pulsars are attributable to their birth in close binary systems. 

The hypothesis of high space velocities of young NSs immediately leads to 
the well-known problem of NS retention in GCs, since the escape velocity 
even in the densest clusters does not exceed several tens of km/s. Most of 
the NSs born through the core collapse of massive stars should have escaped 
from the cluster shortly after its formation. 

In principle, the problem of NSs in GCs can be solved by assuming a non- 
standard mechanism of their formation in GCs, for example, through the 
collapse of white dwarfs in close binary systems induced by mass accretion 
from the secondary component (Helfand et al.1983; van den Heuvel 1984; 
Grindlay and Bailyn 1988; Bailyn and Grindlay 1990). In this scenario, the 
white dwarf must be composed of a mixture of O, Ne, and Mg (Miyaji et 
al.1980); mass accretion onto a carbon-oxygen white dwarf ends with its to- 
tal destruction during the explosive thermonuclear burning that arises when 
its mass approaches the Chandrasekhar limit, producing a type-la super- 
nova (for recent calculations, see, e.g., Dunina- Barkovskaya et al.(2001) and 
references therein). 

Recent evidence suggests that, possibly, not all of the young NSs are born 
with anomalously high space velocities. In particular, this conclusion is 
reached when analyzing evolutionary scenarios that give rise to binary NSs 
like the binary radio pulsar PSR J0737-3039 (Bisnovatyi-Kogan and Tutukov 
2004; Dewi et al.2005; Podsiadlowski et al.2005). In this binary, the mass 
of the young pulsar is 1.25 M Q , the mass of the old (spun up by accretion) 
millisecond pulsar is 1.34 M Q , and the orbital eccentricity is low (e = 0.088). 
This is indicative of a low mass ejected during the explosion of the second 
supernova in this binary and a low kick velocity acquired by the NS dur- 
ing collapse. It was surmised that the NSs formed through the collapse of 
an O-Ne-Mg stellar core triggered by electron captures do not acquire kick 
velocities if these stars are members of binary systems (Podsiadlowski et 
al.2004; van den Heuvel 2004). In all of the remaining cases (a single star 
with M > 8M or a star with a mass > 12 — 14M in a binary), the stellar 
core collapses anisotropically and the forming NSs will acquire significant 
space velocities. As follows from currently available detailed calculations of 
the collapse of O-Ne-Mg stellar cores induced by electron captures (Kitaura 
et al.2005), the mass of the forming NS is 1.25 M . 

Note that the range of initial masses for mainsequence stars in which degener- 
ate O-Ne-Mg cores are formed is not completely clear. In the first calculations 
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(Miyaji et al.1980), this range was 8-12 M , but Iben with coauthors (Ritossa 
et al.1999) subsequently narrowed it to 10-11 M Q ; other evolutionary calcu- 
lations (Pols et al.1998) yield the range 5-8 M Q , depending on the chemical 
composition. Obviously, the particular features of evolutionary calculations 
and, possibly, the rotation of stars are important here. In our calculations, 
we assumed this range to be 8-12 M , since the exact width of this range 
has no effect on the final results due to the presence of a much less certain 
parameter, the initial binary fraction in GCs. 

Simple estimation shows that for the assumed Salpeter initial mass function 
f(M) ~ M~ 2 ' 35 with a minimum mass of 0.1M Q and 100% of binaries in a 
GC with a mass of 5 • 1O 5 M , the expected number of NSs in the cluster, 
according to the hypothesis under consideration, is ~ 1000. This number is 
by two orders of magnitude larger than the expected number of NSs under 
the assumption of a universal kick velocity with a mean value of 400 km/s 
for all NSs, which can be obtained from measurements of the space velocities 
for single radio pulsars (Hobbs et al.2005). 

Clearly, whether the hypothesis about two types of NSs is valid can be es- 
tablished only after elucidating the physical cause of the possible collapse 
anisotropy. However, it seems of interest check whether this hypothesis 
agrees with observational data on the number ofNSs of different types in 
GCs, since the absence of a kick velocity during collapse in some of the NSs 
automatically helps solve the problem of their retention in GCs without in- 
voking other assumptions. This was first pointed out by Pfahl et al.(2002) 
and Podsiadlowski et al.(2004). 

We used an approach to calculating the evolution of compact stars in GCs 
based on the population synthesis method in combination with allowance 
for the dynamical interactions between stars in dense star clusters developed 
at the Sternberg Astronomical Institute of theMoscow State University (see 
Kuranov and Postnov 2004). In this paper, we show that the hypothesis 
under consideration in combination with the universally accepted assumption 
ofNS magnetic field decay at the accretion phase (Bisnovatyi- Kogan and 
Romberg 1974) successfully explains the observed properties of NSs in GCs, 
including their number and distributions in X-ray luminosity (for accreting 
NSs in low-mass X-ray binaries) and spin period (for millisecond pulsars). 



4 



3 THE MODEL 



The fundamental difficulty in describing the evolution of dense star clusters 
lies in the necessity of taking into account the mutual influence of the evo- 
lution of the stellar population on the evolution of the cluster as a whole, 
on the one hand, and the reverse effect of the change in cluster parameters 
on the rate of dynamical interactions between cluster-populating stars, on 
the other hand. The currently available population synthesis codes allow the 
evolution of single and binary stars in the Galaxy to be modeled. There are 
also codes that allow the evolution of self-gravitating systems consisting of a 
large number of star (~ 10 6 ) either of equal masses or in a limited mass range 
to be computed. At present, intensive work is being done to modify and com- 
bine the existing codes into a single code. The international working group 
MODEST (http://www.manybody.org/modest/) was created in recent years 
to solve this problem and substantial progress has been achieved in the above 
field owing to its activity. 

The code presented here implements one of the possible approaches to de- 
scribing and modeling the evolution of binary stars in GCs without resorting 
to direct N-body calculations. The basic method of calculations consists 
in the following. A newly formed binary system (an initial binary or a bi- 
nary formed during the GC evolution as a result of dynamical interactions) 
is placed in a GC where it evolves with allowance made for its interactions 
with the surrounding stars. The evolution of the GC and "background" stars 
is assumed to be known, is identical for each computed binary, and does not 
depend on the evolution of this binary in the GC model used. The evolu- 
tionary track of the binary consists of a finite number of phases. The phases 
change either through the stellar evolution of one of the components (the 
evolution of an "isolated" binary) or through dynamical interactions with 
the surrounding stars. Our code includes the following types of interactions: 

1 . The interactions that lead to changes in the orbital parameters of bi- 
naries (semimajor axis and eccentricity) - flyby; 

2. The close passages that lead to star exchanges; 

3. The passages that lead to the disruption of binary systems; 

4. The interaction between single stars that lead to the formation of binary 
systems (tidal captures). 
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Thus, we assume that the evolution of the cluster and its stellar population 
affects the evolution of the computed binary, but we disregard the reverse 
effect of this binary on the evolution of the cluster and the surrounding stars. 
This approach is somewhat limited (in particularly, it does not allow the 
interactions between binary systems and the interactions with rare objects, 
e.g., black holes, to be taken into account), but it has a number of advantages. 
These include the following: 

f . A fairly high computational speed; 

2. The "clarity" of calculations that allows the degree of influence of par- 
ticular processes on the evolution of binaries to be estimated; 

3. The possibility of calculating an arbitrary (statistically significant) 
number of evolutionary tracks for binary stars within the framework of 
a single GC model. 

3. 1 Initial Distributions and the Evolution of a Binary System in a GC 

An improved version of the scenario machine population synthesis code (Lipunov 
et al.1996) adapted to calculate the evolution in dense stellar systems is used 
to model the detailed evolution of binary systems. We use the following 
standard initial distributions: 

/(log a) = const, max | p^^y[ ^ | < a < 10 3 -R©, (1) 

for the semimajor axes of binary stars and 

/(Mi) oc Mf 2 - 35 , O.1M < Mi < 12OM (2) 

for the mass of the primary component. Here, Rl(MI) is the Roche-lobe 
radius for the primary component. For the component mass ratio, we use a 
power-law distribution: 

f(q) <xq a \ q = M 2 /M 1 < 1. (3) 

Here, a q is the parameter of the distribution in component mass ratio. In 
our calculations presented below, we set a q = 0, i.e., we considered a "flat" 
distribution in component mass ratio. 
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The dipole surface magnetic field B of the forming neutron stars was assumed 
to have a uniform log distribution: 

f{B) dBcxd log B, 10 10 G < B < 10 14 G (4) 

which spans the range of field strengths estimated from observable single 
radio pulsars (except millisecond pulsars) in the absence of magnetic field 
decay. The magnetic field was assumed to decay only for NSs in binary 
systems at the accretion phase. 

In this case, we used an exponential law of magnetic field decay with time 
up to a certain lower limit (van den Heuvel et al.1986): 

(Boe-V*, B>B miQ 

\ B mhl t > r ln(£ /£min) K } 

We took the field decay time scale to be r = 10 7 yr and the minimum 
magnetic field strength to be B min ~ 10 8 G. 

As we noted in the Introduction, NSs acquire a kick velocity during their 
birth. This velocity was modeled as follows: we assumed the direction of the 
kick velocity to be random and isotropically distributed and determined its 
magnitude as 

( = 0, 8M < M < VIM. 

f{v)dv { ocexp(-^)^, M>12M (6) 

for NSs from binaries and as 

f(v)dv ocexpf — - — -J v 2 dv 

for NSs from single stars. In both cases, the distribution parameter w was 
set equal to 250 km/s, in close agreement with the velocity measurements of 
single radio pulsars (Hobbs et al.2005). 

The initial NS spin periods were assumed to be identical for all stars and 
equal to P = 10 ms. 



3.2 Background Stars in GCs 



To properly describe the dynamical interactions of a binary system with other 
GC objects, the code also computes the evolution of single GC stars. The 
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problem is simplified, because information about the distribution of stars only 
in mass (and radius, for tidal captures) at each time will suffice to determine 
the rates of dynamical interaction. All of the remaining parameters of a 
single star at the current time will be required only if it becomes a binary 
component as a result of the exchange interaction or the tidal formation of 
a new binary. This allows us to describe the background of single GC stars 
by creating a discrete set N(m,i, Rj,t), where m and R are, respectively, the 
mass and radius of the star at time t, so that 

i 3 

where N^ otal is the total number of single GC stars at time t. In our calcula- 
tions, we assumed that nm = 45 and nr = 45. The following processes were 
taken into account in the formation of the background of single stars: 

1. The evolution of single stars. We assumed that the evolution of a 
single star was entirely determined by its mass (for a given chemical 
composition): Mj(i) = f(M®,t). Choosing a star from the mass in- 
terval (Mj(t), Mi(t) + Srrii) at time t and using the dependence Mf = 
f~ 1 (Mi(t)), we then randomly draw the initial mass of the star, ac- 
cording to Eq.(2), from the initial interval (Mf,Mf + Sm®). In turn, 
knowledge of the initial mass allows us to calculate the evolution of the 
single star under consideration and to determine all of the remaining 
required parameters at time t. Thus, the set N(rrii, Rj,t) allows the 
stellar population of single GC stars to be completely described. 

2. The space distribution of background stars. 

Mass segregation results in the settling ofmassive stars to the GC cen- 
ter. W e used the Michie-King multimass model (Michie 1963) to de- 
scribe the GC structure. The space density pi of a stellar subsystem 
formed by stars with masses in the mass range (Mj(t), Mj(t) + Srrii(t)) 
(corresponding to the breakdown of the set N(rrii, Rj,t)) is described 
by a power law: 

[ Pc t , r < r c% 

Pi(r) = < P Ci (r/r Ci )- 2 , r Ci < r < r hi (7) 
[ Phi(r/r hi )- 4 , r hi <r< r t 

where r Ci and r/ li are, respectively, the radii of the core and the sphere 
in which half of the mass of the stellar subsystem is contained, and 
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r t is the tidal radius of the cluster. Thus, p^ = p Ci ( r hJ r cJ~ 2 and 
p Ci = M tot ./[8irrl.(r hi — §r c J]. In the case of energy equipartition for 
GC stars, we may assume that 



r, 




Cl 



where m c is the mean mass of the stars in the GC core, r c = y 2?g p 
is the GC core radius, v m is the rms space velocity, and p c is the central 
star number density: 



3. The GC evolution. 

In general, the distribution Pi(r) is a function of time, reflecting the 
GC evolution. In this formulation of the problem, the GC evolution is 
important from the viewpoint of the change in the rate of interaction 
between stars. To take into account the GC evolution (the collapse 
of the cluster core followed by its expansion), we used a power law 
for the changes in parameters r c and v m (t), as follows from numerical 
calculations (see, e.g., Kim et al. 1998). The core collapse at late 
phases is described by the formulas: 



Pc(t) ~ pM (1 - t/Uou)- 1 - 2 ^,^) ~ v m (0) (1 - t/t coll )-°- 12 , (9) 



where t is the current time and t co u is the cluster core collapse time. 

As our numerous calculations show, the core collapse stops as a result 
of the formation of binary stars in the central part of the cluster. The 
succeeding core expansion can also be described by a power law: 



The changes in parameters p c , r c and v m during the cluster evolution are 
schematically shown in Fig.l. We used two GC models with different 
central densities p c and core radii r c at the collapse time (see the table 
and Fig.l). In both models, the collapse time was taken to be t co u = 
7 • 10 9 yr. At the GC formation time, the mass of all cluster stars 
was set equal to 5 • 1O 5 M ; half of all cluster stars were assumed to 
be binary. Roughly speaking, one third and two thirds of the initial 
cluster mass are concentrated in single and binary stars, respectively. 




(8) 



p c (t) ~ t 2 ,V m (t) rsj t 



-0.32 



(10) 
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4. The evaporation of stars from the GC. 

The set iV*(mj, Tj) also changes through the evaporation of background 
stars from the GC and the formation/ disruption of binary systems: 

Nirm^t) = N(m t ,r J ,t-l)-ANl ap (m l ,r 3 )-J2^N 5 b t s (m t ,r 3 ,m k ,r l )), 

k,i 

where N(m,i,rj,t) is the number of GC binaries in the correspond- 
ing phase volume element (Sm,i,Srj), ANl (m,i,rj) is the change in 
N(m,i,rj,t) through the evaporation of stars from the cluster core, 
and the last term describes the decrease in the number of single stars 
through the formation of tidal binaries. In this case, we assumed that 

where the coefficients c evap (m,i) ~ . Lighter stars with lower masses 
mostly go into the halo and their cluster binding energy decreases. This 
behavior is to be expected from the main properties of the slow relax- 
ation via pair collisions. The time scale for this process was estimated 
by Spitzer (1969) (see also Mouri and Taniguchi 2002): 

_ {vl + vlf 12 _ rrn , , 

s 8(67r)V2G 2 m 1 m 2 n 1 ln(A) ^rm^ l} ' 



3.3 Interactions between GC Stars 



The orbital parameters of a binary system moving inside a GC are determined 
by its energy E and angular momentum J. The relative motions of GC 
stars produce continuous fluctuations in the gravitational field; in turn, these 
fluctuations cause the quantities E and J to change. The energy and angular 
momentum can also change through close encounters and if a supernova 
explosion occurred during the evolution of the massive binary component, as 
a result of which the binary instantaneously loses its mass and can acquire an 
additional kick due to collapse anisotropy. When each evolutionary track of 
a binary system is calculated, the changes in E and J with time are tracked; 
for a detailed description, see our previous paper (Kuranov and Postnov 
2004). The rate of interaction between a binary and a subsystem of stars 
i with a number density rij(r) = Pi(r)/m,i at time t is determined by the 
orbit-averaged value 

2 f Va 

Kij{t) = dij x p ^ ^ J m(r, t)v(r, t)dr/v r , 
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where r p and r a are, respectively, the orbital pericenter and apocenter of 
the star in the cluster, v is the space velocity of the binary, v r is the radial 
component, P is the orbital period, 



The subscript j(j = 1. . . 4) corresponds to the type of dynamical interac- 
tion between this binary and single stars and cr^- is the cross section for the 
corresponding process. The probability of a dynamical interaction between 
the binary under consideration and single background stars during time .dyn 
is defined as 



The scenario machine population synthesis method uses a discrete description 
of the evolution of binary components. In a GC, the evolutionary phase 
of a binary changes in time dt = min{dt evo i,Td yn }, i- e. , either through 
the intrinsic evolution of one of the components (the evolution of a binary 
"isolated" from the background) or through its dynamical interactions with 
the surrounding stars. The time dt evo i depends only on the state of the 
binary components at the beginning of each phase and Td yn can be determined 
from Eq. ljll)) . in which the probability of a dynamical interaction p(rd yn ) is 
specified by a random variable uniformly distributed in the interval [0, 1]. 

To choose the type of interaction j and subsystem i of single background 
stars involved in the interaction, let us introduce quantities pij that have 
the meaning of the probability of interaction j with subsystem i of single 
background stars for the binary system under consideration: 



Next, we choose a random number p uniformly distributed in the inter- 
val [0, 1]. We will sumthe elements of matrix p^ (whether over its rows 
or columns) until the sum will exceed the number p. The element i, j in 
which this occurs will define the sought-for values of i and j. This procedure 
ensures an appropriate choice of indices i, j. This procedure is illustrated 
particularly clearly when one or more elements dominate in the matrix p^ 
and when most of the matrix elements are of the same order of magnitude. 






where the total rate of possible interactions is 
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4 RESULTS OF CALCULATIONS 



The results of our calculations for the evolution of binary systems in GCs for 
models A and B are shown in Figs. 2 -8. To obtain statistically significant 
results for each GC model, we computed the evolution of one million binary 
systems by the Monte Carlo method. The background of single stars with 
which the binary system interacts was computed in advance for each GC 
model. The neutron stars that did not escape from the cluster, including 
those formed in binary systems from stars in the mass range 8-12 M & without 
any initial kick velocity, were included in this background. 

In Fig. 2, the total number of binaries with neutron stars is plotted against 
time with a breakdown into subtypes NSs in pairs with main-sequence stars, 
NSs in pairs with white dwarfs, binary NSs, and NS-black hole pairs. The 
binaries retained in the GC are shown in the left panels. The binaries that 
escaped from the GC are shown in the right panels. The dynamical GC core 
collapse time is marked by the vertical dotted line. As might be expected, the 
number of binaries increases after core collapse. The decrease in the number 
of binaries with time for some of the binary types results from component 
mergers (this is particularly clearly seen for the binaries that escaped from 
the cluster) and dynamical disruptions. However, the latter effect causes no 
appreciable decrease in the total number of binary stars inside the cluster, 
but more likely manifests itself as a slowdown in the increase in the number 
of stars with time. The difference in the models of the GC itself is more 
pronounced for the binaries that escaped from the cluster: in the more cen- 
trally condensed cluster B (lower panels), the number of binaries escaping 
from the cluster increases due to more intense dynamical interactions in the 
cluster core. Note also that the number of produced binary NSs proved to 
be very small; to obtain a statistically significant and proper estimate of the 
formation rate of these important astrophysical sources in GCs requires cal- 
culating a larger number of stars and taking into account the collisions with 
background binary stars (Grindlay et al.2005). We plan to perform such 
calculations in the future. 

Figure 3 shows the formation rate of binary systems with neutron stars in 
the case of dynamical interactions between stars inside the GC: exchange 
interactions (an impinging single NS occupies the place of one of the binary 
components) and the formation of binary systems through tidal captures. 
The formation rate through exchange interactions is considerably higher that 
the formation rate through tidal captures for both GC models. In the over- 
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whelming majority of cases, a main-sequence star is the NS companion for 
exchange interactions (the formation rate of binary systems with a white 
dwarf as the secondary component does not exceed lGyr -1 and is not shown 
in the figure). The formation rates of accreting NS observed as low-mass X- 
ray binaries (LMXBs) and millisecond pulsars (MPSRs) formed through the 
same dynamical interactions between stars in the GC (exchange interactions 
and tidal captures) are shown in Fig. 4. 

Figure 5 shows the time evolution of the total number of low-mass X-ray bi- 
naries and millisecond radio pulsars. Different shading indicates the change 
in the number of close binary systems in which accretion onto the NS takes 
place when a main-sequence dwarf and a degenerate low-mass white dwarf 
fill their Roche lobes due to the loss of orbital angular momentum through 
gravitational waves. We see that the number of objects of this type escaped 
from the cluster is comparable to the number of objects left inside the cluster 
and depends weakly on the chosen cluster model. The ratio of the number of 
millisecond radio pulsars to the number of low-mass Xray binaries is shown 
separately in Fig.6. Note that the number of millisecond pulsars spun up 
by accretion in the computed models is approximately half the number of 
accreting NSs. One would think that this situation is inconsistent with the 
observations (recall that one or two strong X-ray sources, many weak X-ray 
sources, and 10-20 millisecond pulsars are observed in a typical GC). How- 
ever , including the initially single NSs at the ejection phase (radio pulsars) 
changes the situation: the ratio of the total number of NSs observed as radio 
pulsars (spun-up pulsars + NSs with initially weak magnetic fields) to the 
number of accreting NS is larger than one (the right panel in Fig.6). It is 
theoretically clear that the initial distributions of NSs in period and magnetic 
field and the evolution of the NS magnetic field are crucially important in 
analyzing these ratios. As regards the observed situation with X-ray sources, 
as we noted above, the situation with weak unidentified X-ray sources re- 
mains uncertain and the number of quiescent LMXBs among the weak X-ray 
sources in GCs can increase with time. 

Figure 7 shows the time evolution of the X-ray luminosity function for GC. 
We see that the luminosity function cannot be described by a single power 
law. The smooth steepening of this function near 10 37 erg s _1 may be due 
to the difference in the regimes of accretion onto the neutron star (Post- 
nov and Kuranov 2005). The mean slope of the X-ray luminosity function, 
d log N/d log L x ~ —1.3, is close to the observed values of the luminosity 
function for low-mass X-ray binaries in galaxies (Fabbiano 2005) before the 
knee at ~ 5 • 10 38 erg s" 1 . In our simplified calculations, we could not obtain 
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stationary X-ray sources with luminosities above the Eddington limit in the 
case of accretion onto the neutron star. 

The period distribution for millisecond pulsars is shown in Fig.8. Recall 
that the initial spin period of all NSs in our calculations was Pq = 10 ms 
(the second peak in the figure corresponds to this value). The presence of 
pulsars with periods shorter than P implies that these stars have passed 
the phase of accretion- driven spin-up in binary systems. We clearly see that 
the fraction of these NSs increases monotonically after the GC core collapse 
stops. The derived peak at millisecond periods is in qualitative agreement 
with the observed period distribution of pulsars in GCs (for the GCs Terzan 
5 and 47 Tuc, see Fig. 3 from Camillo and Rasio 2005). 



5 DISCUSSION AND CONCLUSIONS 

Here, we presented the results of our calculations of the formation and evo- 
lution of neutron stars in dense globular clusters. Our simulations were per- 
formed by the population synthesis of the evolution of binary and single stars 
in the approximation of a persistent background of single stars with which 
they interact dynamically. We considered two GC models with different star 
number densities and different core sizes after dynamical collapse with given 
time evolution of parameters. 

We assumed that the NSs born in binary systems through the collapse of 
degenerate O-Ne-Mg stellar cores induced by electron captures (this was 
thought to be possible at the end of the evolution of stars with initial main- 
sequence masses in the range 8-12 M do not acquire any kick velocity. In 
all of the remaining cases of NS formation (i.e., from single stars with masses 
> 8M or from binary components with masses > 12M ), the newly born 
NSs were assumed to acquire a kick velocity in accordance with a Maxwellian 
distribution with a parameter of 250 km s' 1 . In this scenario, a significant 
fraction of the forming NSs is retained in the cluster. About 1000 NSs re- 
main in a cluster with a mass of 5 • 1O 5 M in which half of its stars were 
initially binary members. In this case, the initial masses of such NSs must 
be 1.25 M Q . Observations of the GC 47 Tuc (Heinke et al.2005) show that 
the space distribution of millisecond pulsars in this cluster corresponds to a 
mean NS mass of ~ 1.39M Q . Since ~ O.1M is required to spin up NSs to 
millisecond periods at the accretion phase (Lipunov and Postnov 1984), these 
measurements can serve as an indication of an initial NS mass of 1.3M . 
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Our modeling of the number ofmillisecond pulsars spun up by accretion in 
binary systems showed good agreement with the observational data: for the 
cluster models under consideration, we obtained ~ 5 — 10 millisecond pulsars 
and approximately the same number of such objects are dynamically ejected 
from the cluster. The period distribution of spun-up pulsars modeled under 
the assumption of accretion-induced NS magnetic field decay to a minimum 
of 10 8 G is also in qualitative agreement with the observational data (Camillo 
and Rasio 2005). 

The situation with accreting NSs is more complex. In our calculations, we 
found their number to be twice that of spun-up millisecond pulsars. The ob- 
served situation is opposite: there are a few lowmass X-ray binaries among 
the identified sources in the GC (Verbunt and Bassa 2004). Note, however, 
that recently the so-called quiescent low-mass X-ray binaries have been in- 
creasingly identified among weak X-ray sources in GCs with luminosities 
L x ~ 10 31 - 10 32 erg s" 1 (see, e.g., Caskett et al. (2005) in Terzan 5 and 
Heinke et al.(2005) in 47 Tuc). The number of modeled binaries of this type 
will probably also decrease when we include the binary-binary interactions, 
which have not yet been included in our model. These types of interaction 
can reduce the binary fraction in GCs with time (for recent calculations, 
see Ivanova et al.2005). Indeed, when interacting with single stars, a binary 
system can be disrupted only if the energy E = of the impinging 

star exceeds the absolute value of the binding energy \E b \ — G " 1 ^" 12 for the 
binary. Here, Voo is the relative velocity of the impinging star (infinitely far 
from the binary); mi,m2, and m are the masses of the binary components 
and the impinging star; and a is the semimajor axis of the binary. Clearly 
, for "rigid" binaries (\E b \ > mV 2 /2, V is the stellar velocity dispersion in 
the GC), the disruptions are more efficient when they interact with binary 
systems. Since we did not consider the interactions between binary systems 
in our calculations, the disruption rate of rigid binaries was definitely under- 
estimated. 

The constructed X-ray luminosity function for accreting NSs in GCs has a 
mean slope, d log N/d log L x ~ —1.3, close to the observed value for lowmass 
X-ray binaries in galactic bulges and elliptical galaxies. This universality 
of the luminosity function for X-ray sources pointed out by Gilfanov (2004) 
can be explained by common properties of the accretion in low-mass X-ray 
binaries (Postnov and Kuranov 2005). 

Note that all of the numerical estimates in our calculations depend signif- 
icantly on the assumed binary fraction in the GC. This parameter may be 
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important for the evolution of the cluster as a whole (Fregeau et al.2005). 
Our calculations do not allow this parameter to be separated from other pa- 
rameters (such as the mass range for stars with degenerate ONe- Mg cores, 
the total number of stars in the GC, etc.). Nevertheless, reasonable agree- 
ment between our calculations and the observational data suggests that our 
main model assumptions about the evolution of stars in the GC are valid. 

Our description of the evolution of binary and single stars in GCs is not 
self-consistent, since it does not take into account some of the dynamical 
processes, such as binary-binary interactions; we also disregarded the reverse 
effect of binaries on the background of single stars by assuming the latter to 
be replenished to the previous level after each interaction. These processes 
are important in calculating the formation of NS+NS pairs inGCs (Grind- 
lay et al. 2005). We plan to take into account these effects in subsequent 
papers. Nevertheless, the main result of our work will not change greatly: 
the hypothesis about the birth of some of the NSs in binary systems without 
any kick velocity allows us to obtain in a natural way a sufficient number 
of neutron stars in globular clusters without invoking an artificial retention 
factor, which becomes unavoidable for explaining the observations if any NS 
during its birth is assumed to acquire a high kick velocity. 
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Table 1. GC model parameters used in our calculations. The central star density 
n c and the core radius r c are given at the core collapse time. 





Initial 


Initial 


Central 


Core 


Collapse 


Model 


GC mass 


binary 


density 


radius 


time, 




M 


fraction 


n c pc~ 3 


r c , pc 


Myr 


A 


5 x 10 5 


0.5 


10 5 


0.3 


7000 


B 


5 x 10 5 


0.5 


10 6 


0.1 


7000 



Figure 1. Changes in cluster parameters: central density p c (pc~ 3 ), core radius r c 
(pc), and core star velocity dispersion v c (km s _1 during the GC evolution. The 
GC core collapse time is indicated by the vertical line (7 Gyr). 

Figure 2. Time evolution of the total number of binaries with neutron stars as 
a function of the type of secondary component: neutron stars in pairs withmain- 
sequence stars (NS + MS), neutron stars with white dwarfs (NS +WD), and binary 
neutron stars (NS + NS). The number of binaries inside the cluster (upper left 
panel) and those escaped from the cluster (upper right panel) for model A and the 
number of binaries insider the cluster (lower left panel) and those escaped from 
the cluster (lower right panel) formodel B are shown. The GC core collapse time 
is indicated by the vertical dotted line (7Gyr). 
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Figure 3. Formation rates of binaries with neutron stars for dynamical interactions 
between stars: the formation rates (Gyr -1 ) in the case of exchange interactions 
between binary and single stars for model A (upper left panel) and model B (upper 
right panel) and the formation rates of binaries through tidal captures for model 
A (lower left panel) and model B (lower right panel) are shown. The GC core 
collapse time is 7 Gyr. 

Figure 4. Formation rates of low-mass X-ray binaries (LMXBs) and millisecond 
pulsars (MPSRs) in the case of interactions between stars: the formation rates 
(Gyr -1 ) in the case of exchange interactions between binary and single stars for 
model A (upper left panel) and model B (upper right panel) and the formation 
rates of binaries through tidal captures for model A (lower left panel) and model 
B (lower right panel) are shown. The GC core collapse time is 7 Gyr. 

Figure 5. Number of accreting neutron stars (as a function of the type of 
Roche-lobe-filling component: (Na + MSGW) in a pair with a main-sequence 
star (LMXBs) and (Na+WD) in a pair with a white dwarf and millisecond pulsars 
(single MPSRs in binaries). The numbers of systems inside the cluster (upper left 
panel) and those escaped from the cluster (upper right panel) for model A and the 
numbers of systems inside the cluster (lower left panel) and those escaped from 
the cluster (lower right panel) for model B are shown. The GC core collapse time 
is indicated by the vertical dotted line (7 Gyr). 

Figure 6. Left panel: ratios of the total number of millisecond pulsars (single + 
binary members) to the number of low-mass X-ray binaries in GC vs. time for 
models A (solid line) and B (dashed line). Right panel: ratios of the total number 
of neutron stars at the ejection phase (single + binary members) to the number of 
low-mass X-ray binaries in GC vs. time for models A (solid line) and B (dashed 
line). 

Figure 7. Fig. 7. Time evolution of the X-ray luminosity function for accreting 
neutron stars in GC formodels A (left panel) and B (right panel), in units of 10 38 
erg s _1 . 

Figure 8. Histograms illustrating the time evolution of the spin period distribution 
for pulsars spun up previously at the accretion phase in GC for models A (left 
panel) and B (right panel). 
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